%Script to validate teh candidate models for IGRF12

%% plot data and secular variation

S2 = dir('/Users/manojnair/data/obs_mag_data/BGS3/*.nc');

% hgdm coefficient files


S = dir('/Users/manojnair/projects/hdgm2015_validation/models/reformatted/*.txt');

for i = 1:length(S),
    S = setfield(S,{i},'dir','/Users/manojnair/projects/hdgm2015_validation/models/reformatted/');
    S = setfield(S,{i},'nmax',15);
    sv = sscanf(S(i).name, 'hdgm_%d');
    S = setfield(S,{i},'svYear',sv);
end;

S(12).svYear = 2014.2;
S(13).svYear = 2014.2;
S(14).svYear = 2014.23;
S(15).svYear = 2014.23;

lengthprevious = length(S);

S1 = dir('/Users/manojnair/projects/igrf12_validation/SV-2015-2020-E.txt');
S = setfield(S,{lengthprevious+1},'dir','/Users/manojnair/projects/igrf12_validation/');
S = setfield(S,{lengthprevious+1},'name',S1(1).name);
S = setfield(S,{lengthprevious+1},'nmax',8);
S = setfield(S,{lengthprevious+1},'svYear',2014.3);

%models


%selected_obs = cell2mat(textread('/Users/manojnair/projects/wmm2015_validation/selected_obs.txt','%s'));

time_array_hour = (datenum(2007,1,1,0,0,30): (1/(24)): datenum(2015,12,31,23,59,30))';

for i = 1:length(S2),
    spline_obs(i,:) = S2(i).name(1:3);
end;

%[data_obs, ia, ib ] = unique(filestring,'rows');

%[spline_obs, ia, ib]= intersect(data_obs, selected_obs, 'rows');

% spline_obs = filestring; % using all observatories

% plot

% f = figure(1);
% set(f,'Position',[40          19        1557         937]);


annual_time_1 = datenum(2010,1:60,15); % 0 year to end year -1
annual_time_2 = datenum(2011,1:60,15); % 1 year to end year
actual_sec_time = datenum(2010,7:66,17); % dates to plot data
annual_time = datenum(2010,1:72,15);

n1 = length(datenum(1995,1,1,0,0,30): (1/(24)): datenum(2015,12,31,23,59,30))';
n2 = length((datenum(1995,1,1,0,0,30): (1/(24)): datenum(2006,12,31,23,59,30))');
n3 = length((datenum(1995,1,1,0,0,30): (1/(24)): datenum(2013,12,31,23,59,30))')-n2;
n4 = length((datenum(1995,1,1,0,0,30): (1/(24)): datenum(2014,12,31,23,59,30))')-n2;

i = 1;
for nobs = 1:length(spline_obs),
    
    ncfname = ['/Users/manojnair/data/obs_mag_data/BGS3/' spline_obs(nobs,:) '_1995_2015.nc'];
    
    
    %display(ncfname);
    
    [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj] = read_geomag_netcdf(ncfname, n2, (n1-n2), 0);
    
    if any(x_data) && sum(~isnan(x_data(n3:n4))) > 8760/2,
        
        
        
        % fit splines
        L = isnan(x_data);
        
        data_array = x_data(~L);
        first_index = find(~L, 1, 'first');
        
        last_index = find(~L, 1, 'last');
        
        time_array = time_array_hour(~L);
        
        b1 = min(time_array):365:max(time_array)+10;
        
        if length(b1) > 1,
            
            sp_x=spline(b1,x_data(~L)'/spline(b1,eye(length(b1)),time_array'));
            sp_y=spline(b1,y_data(~L)'/spline(b1,eye(length(b1)),time_array'));
            sp_z=spline(b1,z_data(~L)'/spline(b1,eye(length(b1)),time_array'));
            
            if  obj.geospatial_lon > 180,
                obj.geospatial_lon = obj.geospatial_lon - 360;
            end;
            
            
            
            geoc_lat = geod2geoc(obj.geospatial_lat, 0, 'WGS84'); % geocentric lat
            geoc_r = ( geocradius(geoc_lat,'WGS84') - 6371.2*1000 )/1000; % geocentric alt in km
            
            % secular variation is linear and same through out the model so
            % no need to calculate it through the loop.
            
            
            geoc_lat = geod2geoc(obj.geospatial_lat, 0, 'WGS84'); % geocentric lat
            geoc_r = ( geocradius(geoc_lat,'WGS84') - 6371.2*1000 )/1000; % geocentric alt in km
            
            obs_data(i).name = spline_obs(i,:);
            
            
            for kk = 1:length(S),
                
            obs_data(i).model_sv(kk).data = magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2015.0,[S(kk).dir S(kk).name],S(kk).nmax);

            end;
            
            obs_data(i).dat_2014 = sum(~isnan(x_data(n3:n4)));
           
            for kk = 1:length(S),
            
            obs_data(i).sec_x_obs(kk) =  ppval(sp_x,datenum(S(kk).svYear + 0.5 ,0,0)) - ppval(sp_x,datenum(S(kk).svYear - 0.5 ,0,0)) ;
            obs_data(i).sec_y_obs(kk) =  ppval(sp_y,datenum(S(kk).svYear + 0.5 ,0,0)) - ppval(sp_y,datenum(S(kk).svYear - 0.5 ,0,0)) ;
            obs_data(i).sec_z_obs(kk) =  ppval(sp_z,datenum(S(kk).svYear + 0.5 ,0,0)) - ppval(sp_z,datenum(S(kk).svYear - 0.5 ,0,0)) ;
            
            end;
            
            fprintf('%s \n', spline_obs(nobs,:));
            
%             subplot(311);
%             plot(time_array_hour, x_data);
%             datetick('x','yyyy','keeplimits')
%             subplot(312);
%             plot(time_array_hour, y_data);
%             datetick('x','yyyy','keeplimits')
%             
%             subplot(313);
%             plot(time_array_hour, z_data);
%             datetick('x','yyyy','keeplimits')
            
            obs_data(i).lat = obj.geospatial_lat;
            obs_data(i).lon = obj.geospatial_lon;
            
            title(spline_obs(i,:));
         
            
            clear  sp_x sp_y sp_z;
            %pause;
            clf;
            
            i = i+1;
        end;
        
    end;
    
end;


nobs = i-1;


%% print ms results


for i = 1:nobs,
    
    for j = 1:length(S),
        
        model_sv_x(i,j) = obs_data(i).model_sv(j).data(1);
        model_sv_y(i,j) = obs_data(i).model_sv(j).data(2);
        model_sv_z(i,j) = obs_data(i).model_sv(j).data(3);
        
        obs_sv_x(i,j) = obs_data(i).sec_x_obs(j);
        obs_sv_y(i,j) = obs_data(i).sec_y_obs(j);
        obs_sv_z(i,j) = obs_data(i).sec_z_obs(j);
    end;
end;
    
for i = 1:length(S)
    
fprintf('%s %5.2f %5.2f %5.2f\n', S(i).name,rms((model_sv_x(:,i) - ...
    obs_sv_x(:,i))-mean(model_sv_x(:,i) - obs_sv_x(:,i))), ...
    rms((model_sv_y(:,i) - obs_sv_y(:,i))-mean(model_sv_y(:,i) - obs_sv_y(:,i))), ...
    rms((model_sv_z(:,i) - obs_sv_z(:,i)) - mean(model_sv_z(:,i) - obs_sv_z(:,i))));
end;


for i = 1:length(S)
    
fprintf('%s %5.2f %5.2f %5.2f\n', S(i).name, mean(model_sv_x(:,i) - obs_sv_x(:,i)), ...
    mean(model_sv_y(:,i) - obs_sv_y(:,i)), ...
    mean(model_sv_z(:,i) - obs_sv_z(:,i)));
end;
  



%% print mean results


fprintf('A %5.2f %5.2f %5.2f\n', mean(sec_A_x - [obs_data.sec_x_obs_A]), ...
    mean(sec_A_y - [obs_data.sec_y_obs_A]), mean(sec_A_z- [obs_data.sec_z_obs_A]));

fprintf('B %5.2f %5.2f %5.2f\n', mean(sec_B_x - [obs_data.sec_x_obs_B]), ...
    mean(sec_B_y - [obs_data.sec_y_obs_B]), mean(sec_B_z- [obs_data.sec_z_obs_B]));

fprintf('E %5.2f %5.2f %5.2f\n', mean(sec_E_x - [obs_data.sec_x_obs_E]), ...
    mean(sec_E_y - [obs_data.sec_y_obs_E]), mean(sec_E_z- [obs_data.sec_z_obs_E]));

